In [ ]:
df_both = df_both.loc[(df_both["iHep"] != 0.0)& (df_both["HepG2"] != 0.0)]
clustermap = sns.clustermap(df_both.fillna(0), figsize =(10,20))
import matplotlib.pyplot as plt
import cobra
from cobra.io import validate_sbml_model
import importlib
from xml.etree import ElementTree
import utils.Model_correction as mc
import sys
import utils.model_maj as mj
import utils.viz_utils as vu
#import cplex
import cbmpy
import pandas as pd
import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
from tqdm import tqdm
from itertools import cycle
import seaborn as sns
importlib.reload(mj)
importlib.reload(vu)
importlib.reload(mc)
HepG2, errors = validate_sbml_model("../models_storage/Hep_G2_v7.xml")
r = HepG2.reactions.HMR_0030
r.compartments.add("C_x")
r.boundary
iHep, errors = validate_sbml_model("../models_storage/iHep_v3.xml")
i = 0
for S_r in iHep.boundary :
if S_r.bounds == (-1000.0,1000.0) :
i+=1
else :
vu.print_reactions(S_r)
print(S_r.bounds, "\n\n")
print(f"{i}/{len(iHep.boundary)}")
iHep_ex = mj.maj("iHep_v2.xml", "Hep_G2_v6.xml", "../models_storage/", bounds_check = False, genes_id_copy= False, alt_gene_ids= False, metab_id_check= False, bounds_value_check= False, subsystem_copy= False,add_exchanges=False, get_names=True)
iHep_ex.reactions.HMR_4281.name
cobra.io.write_sbml_model(HepG2, "../models_storage/Hep_G2_v7_mod.xml")
# Bidouillage pour ajouter des subsystems aux réactions d'échange. (Va ajouter "exchange reactions" normalement.)
for r in HepG2.reactions :
if "EX_" in r.id :
m = [m_id for m_id, _ in HepG2.reactions.get_by_id(r.id).metabolites.items()]
sub = [r.subsystem for r in m[0].reactions if "EX_" not in r.id]
r.subsystem = sub[0]
HepG2_C = HepG2.copy()
HepG2.reactions.EX_m01965x.bounds = (-565.0,-565.0) #Glucose exchange, set to the maximal input value observed with FVA
#iHep.reactions.EX_m01965x.bounds = (-1000.0, -1000.0)
iHep.reactions.HMR_4281.bounds = (0.0,0.0)
HepG2.reactions.HMR_3883.bounds = (0.0,0.0) #Serine --> Pyruvate exchange
HepG2.reactions.EX_m02630x.bounds = (0.0,0.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
HepG2.reactions.HMR_4316.bounds = (-1000.0,1000.0) # Glucose --> D-Glucitol
HepG2.reactions.EX_m02896x.bounds = (0.0,1000.0) # Serine intake
HepG2.reactions.EX_m01682x.bounds = (-1000.0,0.0) # Glucitol secretion blocked. Kept the intake just in case.
HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.HMR_4930.bounds = (-1000.0,1000.0) # Pyruvate transfer from cytoplasm to peroxysome
#HepG2.reactions.HMR_4281.bounds = (0.0,0.0) # Fermentation in peroxysome
# ADI-PEG20 treatment
HepG2.reactions.EX_m01365x.bounds = (0.0,0.0) # Arginine intake
HepG2.reactions.HMR_3813.bounds = (0.0,0.0) # Arginosuccinate synthase
#HepG2.reactions.EX_m02403x.bounds = (0.0,1000.0) #Lactate exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
#HepG2.reactions.EX_m01840x.bounds = (0.0,0.0) #Fructose exchange
HepG2.reactions.HMR_4281.bounds = (-1000.0,1000.0) # Pyruvate fermentation in peroxysome
HepG2.reactions.HMR_4281.name = "Lactic Fermentation"
HepG2.reactions.EX_m02630x.bounds = (0.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02630x.bounds = (-1000.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
#HepG2.reactions.EX_m01965x.bounds = (-1000,1000.0) #Glucose exchange
#HepG2.reactions.EX_m01286x.bounds = (0.0,0.0) #ADP-Glucose exchange
#HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.EX_m01743x.bounds = (-1000.0,1000.0) #D-Ribulose exchange
#HepG2.reactions.EX_m02453x.bounds = (0.0,0.0) #Mannose exchange
#HepG2.reactions.EX_m02470x.bounds = (0.0,0.0) #Methanol exchange
# --> After Methanol cut, Glucose uptake flux went from 0.0 to -190.0
# Without O2, glucose uptake went down to -23
iHep.objective = "biomass_components"
sol = iHep.optimize()
sol
HepG2.objective = "biomass_components"
val = HepG2.optimize()
val.objective_value
HepG2_C.objective = "biomass_components"
val = HepG2_C.optimize()
val
[(m.name, m.id) for m in HepG2.reactions.biomass_components.metabolites]
importlib.reload(vu)
vu.print_exchanges(HepG2, filter = "all")
vu.run_parcours(HepG2.reactions.EX_m01965x, HepG2)
#m01682c glucitol
for reaction, flux in f.items() :
vu.print_reactions(HepG2.reactions.get_by_id(reaction), flux)
print(f"FLUX : {flux} --- ID : {HepG2.reactions.get_by_id(reaction).id} --- compartment : {HepG2.reactions.get_by_id(reaction).compartments} \n\n---\n\n")
vu.print_reactions(HepG2.reactions.biomass_components)
HepG2.reactions.EX_m02630x.bounds = (-1000.0,1000.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
#HepG2.reactions.EX_m01965x.bounds = (-1000,1000.0) #Glucose exchange
#HepG2.reactions.EX_m01286x.bounds = (0.0,0.0) #ADP-Glucose exchange
#HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
#HepG2.reactions.EX_m01743x.bounds = (-1000.0,1000.0) #D-Ribulose exchange
#HepG2.reactions.EX_m02453x.bounds = (0.0,0.0) #Mannose exchange
#HepG2.reactions.EX_m02470x.bounds = (0.0,0.0) #Methanol exchange
# --> After Methanol cut, Glucose uptake flux went from 0.0 to -190.0
# Without O2, glucose uptake went down to -23
from cobra.flux_analysis import flux_variability_analysis
fva = flux_variability_analysis(HepG2)
HepG2.objective = "biomass_components"
indispensable = []
unisens = []
reversible = []
unknown = []
other = []
for reaction_id, flux_values in fva.iterrows() :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
difference = abs(max_flux - min_flux)
if max_flux > 0.0 and min_flux > 0.0 :
f = max_flux
elif max_flux > 0.0 and min_flux < 0.0 :
f = 0.0
elif max_flux < 0.0 and min_flux < 0.0 :
f = min_flux
else :
f = 0.0
if difference < 1e-9 and min_flux + max_flux != 0.0 and abs(min_flux + max_flux) > 0.01 :
indispensable.append((reaction_id, flux_values, f))
elif min_flux < 0.0 and max_flux > 0.0 :
reversible.append((reaction_id, flux_values, f))
elif min_flux > 0.0 and max_flux > 0.0 or min_flux < 0.0 and max_flux < 0.0 :
if difference > 1e-9 :
unisens.append((reaction_id, flux_values, f))
elif min_flux == max_flux == 0.0 :
unknown.append((reaction_id, flux_values, f))
else :
other.append((reaction_id, flux_values, f))
print("\n--------INDISPENSABLE--------\n")
for reaction_id, flux_values, f in indispensable :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX")
print(f"\nID : {reaction_id} -- NAME : {HepG2.reactions.get_by_id(reaction_id).name}\n\n")
print("\n--------UNISENS--------\n")
for reaction_id, flux_values, f in unisens :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX")
print(f"\nID : {reaction_id} -- NAME : {HepG2.reactions.get_by_id(reaction_id).name}\n\n")
print("\n--------REVERSIBLE--------\n")
for reaction_id, flux_values, f in reversible :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX")
print(f"\nID : {reaction_id} -- NAME : {HepG2.reactions.get_by_id(reaction_id).name}\n\n")
if len(unknown) >0 :
print("\n--------UNKNOWN--------\n")
for reaction_id, flux_values, f in unknown :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX")
print(f"\nID : {reaction_id} -- NAME : {HepG2.reactions.get_by_id(reaction_id).name}\n\n")
if len(other) >0 :
print("\n--------OTHER--------\n")
for reaction_id, flux_values, f in other :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
vu.print_reactions(HepG2.reactions.get_by_id(reaction_id), f)
print(f"\nMIN FLUX -- {min_flux} === {max_flux} -- MAX FLUX")
print(f"\nID : {reaction_id} -- NAME : {HepG2.reactions.get_by_id(reaction_id).name}\n\n")
def cut_and_optimize(model:cobra.core.model.Model, reaction_id:str, objective:str = "biomass_components") :
reaction_to_cut = model.reactions.get_by_id(reaction_id)
original_bounds = reaction_to_cut.bounds
reaction_to_cut.bounds = (0.0,0.0)
model.objective = objective
solution = model.optimize()
objective_value = solution.objective_value
if objective_value:
return objective_value
else :
reaction_to_cut.bounds = original_bounds # Backtrack if no solution is found
values = []
reactions_cut = []
for reaction_id, flux_values in fva.iterrows() :
min_flux = flux_values["minimum"]
max_flux = flux_values["maximum"]
if max_flux == 0.0 and min_flux == 0.0 :
pass
elif max_flux == 0.0 or min_flux == 0.0 :
val = cut_and_optimize(HepG2, reaction_id)
if val == None :
values.append(0.0)
else :
values.append(val)
reactions_cut.append(reaction_id)
else :
pass
for value, id in zip(values, x) :
if value == 0.0 and id < 1200 :
print(HepG2.reactions.get_by_id(reactions_cut[id]).name)
x = [i for i in range(len(values))]
plt.rcParams["figure.figsize"] = (20,10)
fig, ax = plt.subplots()
ax.plot(x, values)
fig.show()
unknown_ids = [i[0] for i in unknown]
for id in unknown_ids :
r_to_change = HepG2.reactions.get_by_id(id)
r_to_change.bounds=(0.0,0.0)
print(f"\nchanged {id}'s bounds")
with open("zero_reactions.txt", "w") as out :
for id in unknown_ids :
out.write(f"{id}\n")
for r in HepG2.reactions :
if len(r.name) > 1 :
print(r.name)
df_Ci, df_Ri, df_Si, df_Mi, df_Pi, df_Xi, df_Li, df_Gi, df_Ni = vu.build_reaction_df(iHep)
df_C, df_R, df_S, df_M, df_P, df_X, df_L, df_G, df_N = vu.build_reaction_df(HepG2)
subS_HepG2 = vu.get_subsystem_fluxes(vu.build_reaction_df(HepG2))
subS_iHep = vu.get_subsystem_fluxes(vu.build_reaction_df(iHep))
subS_G2C = vu.get_subsystem_fluxes(vu.build_reaction_df(HepG2_C))
df_both = pd.concat([subS_HepG2, subS_G2C, subS_iHep], axis = 1)
df_both.columns = ["HepG2","HepG2_copy", "iHep"]
df_both = df_both.loc[(df_both["iHep"] != 0.0)& (df_both["HepG2"] != 0.0)]
clustermap = sns.clustermap(df_both.fillna(0), figsize =(10,20))
Heatmap : Représentation sous forme d'une heatmap des flux de chaque sous-systèmes actif dans trois modèles de cellule du foie différentes :
On s'attend à observer une distribution des flux dirigée vers la production d'énergie, de nucléotides, d'acides aminés ou d'acides gras pour les modèles cancéreux.
La glycolyse est effectivement plus active dans le modèle HepG2 corrigé. \ On observe également une activité importante du métabolisme du folate dans le modèle HepG2 (métabolisme des nucléotides)\ Cependant, l'activité du métabolisme des pentoses phosphates reste bien trop faible, il en va de même pour le métabolisme du pyruvate.
On s'attend à observer les différences suivantes dans le modèle cancéreux par rapport au modèle sain :
--> Résultats :
barplots_dict = vu.compartment_fluxes_barplots(iHep, HepG2)
title = "Cytoplasm reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_C,HepG2,title).show(renderer = 'notebook')
title = "Cytoplasm reactions -- <b>HEALTHY MODEL <br />"
vu.plot_treemap(df_Ci, iHep, title).show(renderer = "notebook")
barplots_dict["C_c"][0]
title="Mitochondrial reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(df_Mi, iHep, title).show(renderer = "notebook")
title="Mitochondrial reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_M, HepG2, title).show(renderer = "notebook")
barplots_dict["C_m"][0]
title="Peroxysomal reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(df_Pi, iHep, title).show(renderer = "notebook")
title_text="Peroxysomal reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_P, HepG2, title).show(renderer = "notebook")
barplots_dict["C_p"][0]
def get_exchanges_reactions(target_model) :
Ex_reactions_to_add = []
target_model.compartments["C_x"] = "C_x"
for C_s_reaction in tqdm(target_model.boundary) :
if len(C_s_reaction.products) == 0 :
if len(C_s_reaction.reactants) == 1 :
# Création du métabolite externe
x_metabolite_id = C_s_reaction.reactants[0].id[:-1] + "x" # Creating external metabolite ID
x_metabolite = C_s_reaction.reactants[0].copy()
x_metabolite.id = x_metabolite_id
target_model.add_metabolites([ \
cobra.core.metabolite.Metabolite( \
x_metabolite_id, \
name=x_metabolite_id, \
compartment='C_x')])
# Création de la réaction externe.
C_x_reaction_id = C_s_reaction.id = "EX_" + x_metabolite_id # Copying C_s reaction to make it an EX_ reaction
# Ajout du métabolite externe à la réaction S
C_s_reaction.add_metabolites({x_metabolite : 1.0})
# Ajout du métabolite externe à la réaction X, et enlèvement du métabolite S.
try :
C_x_reaction = target_model.add_boundary(x_metabolite, type="exchange",reaction_id=C_x_reaction_id, lb=0.0, ub=1000.0,sbo_term="SBO:0000627")
Ex_reactions_to_add.append(C_x_reaction)
except ValueError :
continue
elif len(C_s_reaction.reactants) == 0 :
if len(C_s_reaction.products) == 1 :
# Création du métabolite externe
x_metabolite_id = C_s_reaction.products[0].id[:-1] + "x"
x_metabolite = C_s_reaction.products[0].copy()
x_metabolite.id = x_metabolite_id
target_model.add_metabolites([ \
cobra.core.metabolite.Metabolite( \
x_metabolite_id, \
name=x_metabolite_id, \
compartment='C_x')])
# Création de la réaction externe.
C_x_reaction_id = C_s_reaction.id = "EX_" + x_metabolite_id
C_s_reaction.add_metabolites({x_metabolite : -1.0})
try:
C_x_reaction = target_model.add_boundary(x_metabolite, type="exchange",reaction_id=C_x_reaction_id, lb=-1000.0, ub=0.0,sbo_term="SBO:0000627")
Ex_reactions_to_add.append(C_x_reaction)
except ValueError :
continue
else :
pass
target_model.add_reactions(Ex_reactions_to_add)
iHep_to_fix, errors = cobra.io.validate_sbml_model("../models_storage/iHepatocytes2322_Cobra.xml")
help(cobra.core.reaction.Reaction)
438 boundary reactions originally
boundary_reactions = iHep_to_fix.boundary
reactions_to_add = []
for b_r in boundary_reactions :
for metab, sto in b_r.metabolites.items() :
# Creating boundary metabolite
lb = b_r.lower_bound
ub = b_r.upper_bound
new_metab_x = cobra.core.metabolite.Metabolite(metab.id[:-1]+"x", name = metab.id[:-1]+"x", compartment = "C_x")
if sto > 0.0 :
b_r.add_metabolites({new_metab_x: -1.0}) # Intake
else :
b_r.add_metabolites({new_metab_x : 1.0}) # Secretion
new_reaction_x = cobra.core.reaction.Reaction(id = "EX_" + new_metab_x.id, name = "EX_" + new_metab_x.id, subsystem= "Exchange reactions", lower_bound= lb, upper_bound = ub)
new_reaction_x.add_metabolites({new_metab_x : 1.0})
reactions_to_add.append(new_reaction_x)
iHep_to_fix.add_reactions(reactions_to_add)
cobra.io.write_sbml_model(iHep_to_fix, "../models_storage/iHep_v2.xml")
len(HepG2.boundary)
list(rtest.genes)[0].id